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ABSTRACT 

Inherent instabilities in the radiative transfer for astrophysical masers have been recognized and calculated 
in the linear maser idealization in our previous investigations. The same instabilities are now shown to occur 
in the more realistic, three-dimensional geometries. Fluctuations in the emergent flux result and may be 
related to the observed fluctuations in the radiative flux from the 1665 MHz OH masers that have been 
reported to occur on timescales as short as 1000 s. The time-dependent differential equations of radiative 
transfer are solved numerically for three-dimensional astrophysical masers. Computations are performed for 
spherical and elongated (rectangular parallelepiped) geometries. 

Subject headings: instabilities — masers — radiative transfer 

1. INTRODUCTION 2. BASIC METHODS 

Clegg & Cordes (1991) have reported fluctuations of ~ 10% 
in the flux of radiation from interstellar 1665 MHz OH masers 
on timescales down to 1000 s. Calculations have been per- 
formed in which inherent instabilities in the radiative transfer 
are found to be a likely cause for the observed fluctuations 
(Scappaticci & Watson 1992a, b; hereafter SW1 and SW2). The 
occurrence of the instabilities has, however, been demonstrated 
only in the usual idealization of a linear maser. In the linear 
maser idealization, all rays of maser radiation follow the same 
path through the maser. In actual astrophysical masers, the 
rays are separated by distances of 10 13 - 10 14 cm— the observed 
dimensions of the masers. For significant fluctuations in the 
observed maser flux, the entire cross section of the masing 
region (or at least a large fraction) must respond collectively to 
the instabilities. It may not be evident from the investigations 
to date that such collective behavior is a consequence of the 
instabilities that have been calculated for linear masers. The 
purpose of this paper is the demonstration of these instabilities 
for the radiative transfer in three-dimensional masers as rep- 
resented by spheres and cylinders. Although the range of the 
values for the relevant parameters in which the masers are 
unstable (both linear and three-dimensional) has so far been 
found to be limited to those that are plausible for the 1665 
MHz and other 18 cm OH masers, the basic instability may 
appear in other forms. Rapid time variations of the flux is a 
pervasive feature of astrophysical masers (e.g.. Argon et al. 

1994). 

In SW1 and SW2, the instabilities were established both by a 
linear stability analysis and by numerically integrating the full, 
partial differential equations in space and time. The insta- 
bilities were evident in both for exactly the same range of 
values for the relevant parameters. This indicates that either 
procedure is reliable. For the spheres and cylinders considered 
in this paper, only the results of the numerical integration of 
the full equations are used to examine the instabilities. They 
are found to occur at essentially the same range of parameter 
values as for linear masers — a result that tends to be reassuring 
about the reliability of the calculations presented in § 3. 


The time-dependent equations of radiative transfer in three 
dimensions are solved by generalizing our previous methods 
for a linear maser in a straightforward manner The masing 
transition is treated in the common, two-level approximation 
in which the interaction with other molecular states is approx- 
imated by “phenomenological” pumping rates A and decay 
rates T. With the usual Rayleigh-Jeans approximation for the 
maser radiation, the time-dependent rate equation for the dif- 
ference per unit volume between the populations of the upper 
(u) and lower (l) states of the masing transitions is given by 


1-1 +P' 


p(i, 0 = r(n u - w,)/AA (2) 

is the normalized difference between the populations per mag- 
netic substate of the upper (nj and lower (n,) states at location 
x and time r, AA = A u — A,, A is the Einstein coefficient for 
spontaneous emission, 

«' 1 kT ° A ni 


v is the frequency of the maser transition, and 

P s = T(n u + n,)/AA = (A u + A,)/AA = constant . (4) 

The integral of the normalized brightness temperature f (i, f) 
of the maser radiation over all angles at location jc in equation 
(1) gives the effect of stimulated transitions on the populations. 
Normalized time and spatial coordinates l — t/(Lfc) and x = 
x/L are expressed in terms of the characteristic distance L 
which is the length for cylinders or the diameter for spheres. 
Incident continuum radiation has a brightness temperature T 0 . 
The equation of radiative transfer for the intensity of maser 
radiation expressed as a normalized brightness temperature 
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T (i, f) is then 

dT df (a, hv hv \ 

ds + dl ~ + 2kT 0 P + 2kT 0 Ps ) ’ ( ) 

where s is the normalized distance along the path of the ray of 
radiation and 


t — 



A_ AA 
Av r L 


( 6 ) 


and T (i, T) is the actual brightness temperature divided by T 0 . 

As in the previous calculations, a rectangular distribution of 
width Av in velocity is adopted for the molecular populations 
and an analogous rectangular profile is adopted for the maser 
radiation. Rapid relaxation of the molecular velocities is thus 
assumed. Computations are then performed at only a single 
velocity in equation (1) and frequency in equation (5). These 
approximations were also made in the previous discussions 
and were discussed at that time. Subsequently, computations 
have been performed for a linear maser in which the Maxwel- 
lian distribution for the particle velocities is retained. The dis- 
tribution in frequency for T is then obtained as a result of 
solving equation (5) at a sufficient number of frequencies to 
define the spectral line profile. The effect of these changes on 
the stability of the masers is negligible (Emmering & Watson 
1994). 

To perform the computations, we proceed in the following 
manner. Locations are specified at which the population differ- 
ences p(x y t) are to be computed from equation (1). A separate 
group of rays is then specified for each such location in order 
to compute the integral of the intensity over angles that is 
required in equation (1). A group of such rays is thus “ centered 
upon” the location at which the group is used for computing 
p(x y t). That is, each member of the group passes through that 
location and the directions of the rays are chosen for optimal 
efficiency in computing the integral. These directions are, in 
general, different for each location. Equation (5) is used to 
compute the intensities along these rays as a function of time 
and as a function of the distance along the rays by taking steps 
As sb cAt — exactly as was done in SW2. Here there are a large 
number of rays, whereas in SW2 there were only two rays — the 
“forward” and “backward” rays in the linear maser. In all 
cases, we verify that the results are unaltered when changes are 
made in the number of rays, the number of locations at which 
p(x y f) is calculated, and the step size used in solving equation 
(5). For a specified amount of computer time, we find that 
optimal accuracy is achieved when the separation between the 
locations at which the p(x y t) are computed is much greater (at 
least a factor of 10) than the step size used in solving equation 
(5) for the intensities. At intermediate locations, the p(x, t) are 
obtained by interpolation. In the computations for the rec- 
tangular parallelepiped, the p(x y f) are computed at the grid 
points of an 8 x 8 x 16 lattice. For the sphere, the only 
unstable modes that are identified are spherically symmetric. 
The results presented here are thus based on computations in 
which p( x, t) is found only at a number (24) of radial locations. 
Exploratory computations have also been performed in which 
radial symmetry is not imposed. Unstable, nonradial modes 
have not been found. The sphere is inscribed within a 
24 x 24 x 24 cubic grid of lattice points at which the popu- 
lations p( Jf, t) are computed. A sphere having a diameter equal 
to the length of an edge of the cubic grid is created by setting 
the populations p(x y r) and p s to zero outside of this spherical 


volume. As a guideline, the approximate number of rays in a 
group that is considered to perform the angular integral at 
each location in equation (1) is the number of lattice points on 
the surface of the rectangular parallelepiped or on the surface 
of the cube. That is, approximately 640 (for the 8 x 8 x 16 
grid) and 3456 (for the 24 x 24 x 24 grid), respectively. Sym- 
metries can often be utilized to reduce this number. As in SW2, 
the evolution of the system is initiated by instantaneously 
“turning on” the pumping. That is, equations (1) and (5) are 
solved beginning with p( x, l — 0) = 0 and T (x, f — 0) = 1 
everywhere. How the time evolution is begun was found to be 
unimportant in SW2 for identifying instabilities and this start- 
ing condition is computationally most efficient. After a long 
time and regardless of the starting conditions, it was found that 
the radiation from the unstable masers oscillates and that from 
the stable masers becomes constant. Based on this criterion, 
complete agreement is demonstrated in SW2 between the 
results obtained from the time evolution and from the indepen- 
dent, linear stability analysis with regard to which masers are 
stable and which are unstable. 

3. RESULTS 

Computations are performed for spherical masers and for 
elongated masers as represented by rectangular parallelepipeds 
with square cross sections. The ratio of the lengths to the 
widths of these rectangular parallelepipeds is 100. These are 
intended to represent the extremes of the likely geometries. The 
computations are quite time-consuming, and considerable care 
must be taken to assure convergence. Only limited results are 
presented (a) to establish that the basic instability does occur 
for three-dimensional geometries and ( b ) to indicate that the 
values for the parameters at which the instability occurs are 
similar to those at which it occurs in the linear masers. 

The occurrence of permanent oscillations for the emergent 
maser radiation in Figures 1, 2, and 3 establishes that both 



Fig. 1. — The emergent intensity (as represented by the natural log of the 
normalized brightness temperature) vs. time for the ray of radiation that is 
along the axis of an elongated, three-dimensional maser in the form of a 
rectangular parallelepiped with a square cross section. Time is expressed in 
units of Lfc y where L is length of the maser. The maser has a length to width 
ratio of 100. Results are shown for three choices of the relevant parameters. 
Curve b represents an unstable maser; curve c, a stable maser and curve a, a 
barely stable maser. 
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Fig. 2. — The emergent intensity of maser radiation as a function of time, 
with the same meaning as Fig. I, except that the influence of spontaneous 
emission is examined. The importance of spontaneous emission depends upon 
the parameters p s and T 0 . 

elongated and spherical masers are subject to the instability 
identified in SW1 and SW2. The intensities are presented in the 
figures for only a single, representative ray of radiation. For 
other rays, the variations are similar and have the same phase. 
As for linear masers, the instability occurs when TL/c 2 * 1 
(where L is the length of the elongated masers or the diameter 
of the spheres) and causes oscillations with a period of approx- 
imately L/c. Comparison of curves b and c in Figure 1 indi- 
cates that the instability also disappears at other TL/c. The 
parameter /?' here is analogous to P in SW1 and SW2. Just as 
for the linear masers, the instability disappears for larger 
values of this parameter. For the OH masers at 1665 MHz, 
P' = 1.4 x 1(T 5 (T o /100 K) (0.001 s’VH. Plausible values for 
T 0 and T (e.g., SW1 and SW2) are such that P' can be less than 
10“ 5 as is indicated for instability in Figure 1. Spontaneous 
emission was previously found to reduce the tendency for 
instability in linear masers. It has the same effect in three- 
dimensional masers as can be seen by comparing curve b of 
Figure 2 with curve b in Figure 1. The relative importance of 
spontaneous emission decreases with decreasing fractional 
inversion in the pumping — a 1% difference between the 



Fig. 3. — Same as Figs. I and 2, except that here the maser is in the form of a 
sphere for which the characteristic length L is the diameter of the sphere. The 
intensity is shown for a ray that emerges along a diameter. 


pumping rates of the upper and lower states corresponds to 
p t = 50 and is a typical estimate for the 1665 MHz masers as is 
T 0 ~ 10-100 (e.g., SW1 and SW2). The diagonal ray between 
opposite corners of rectangular parallelepipeds is longer, and 
thus tends to be more intense, than the ray along the axis. For 
the aspect ratio (100) of the parallelepipeds in Figures 1 and 2, 
the difference is negligible. This ray is, of course, included in the 
computations. In exploratory computations for cubic and 
other rectangular parallelepipeds with lower aspect ratios, the 
instability persists and the time variations are similar in char- 
acter to those presented in Figures 1 and 2. Obtaining con- 
vergence is more difficult for spherical massers than for 
elongated masers. The two computations are thus presented in 
Figure 3 only to establish the principle that spherical masers 
are also subject to the same instability. The stable solution is 
presented to demonstrate that our computational methods 
also lead to stable spherical masers and that this occurs for 
parameter values at which it would be expected by analogy 
with linear masers (here, at larger /?'). Spherical masers tend to 
be more stable than elongated masers. For the instability of 
spherical masers, and based upon the information from our 
computations, P ' must be smaller by at least a factor of 10 than 
for the elongated masers. 

Application of the calculated instabilities to the observations 
of Clegg & Cordes (1991) has been discussed in detail in SW2. 
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